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Conventional ordering transitions, described by the Landau paradigm, are characterized by the symmetries 
broken at the critical point. Within the constrained manifold occurring at low temperatures in certain frustrated 
systems, unconventional transitions are possible that defy this type of description. While the critical point exists 
only in the limit where defects in the constraint are vanishingly rare, unconventional criticality can be observed 
throughout a broad region of the phase diagram. This work presents a formalism for incorporating the effects of 
such defects within the framework of scaling theory and the renormalization group, leading to universal results 
for the critical behavior. The theory is applied to two transitions occurring within a model of spin ice, and the 
results are confirmed using Monte Carlo simulations. Relevance to experiments, particularly in the spin-ice 
compounds, is discussed, along with implications for simulations of related transitions, such as the cubic dimer 
model and the 0(3) sigma model with "hedgehog" suppression. 
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I. INTRODUCTION 

The conventional understanding of phase transitions is 
based on the Landau paradigm, in which phases are classified 
by their broken symmetries. 1 Central to this approach is the 
order parameter, a local quantity that transforms nontrivially 
under the symmetries, and hence takes a nonzero expectation 
value only in an "ordered", i.e., symmetry-broken, phase. The 
critical behavior near a continuous phase transition, including 
scaling and universality, can be understood, via the renormal- 
ization groupJ33 in terms of the long-wavelength fluctuations 
of the appropriate order parameter. 

The resounding success of this approach, for both thermal 
and quantum phase transitions,^ has caused considerable at- 
tention to be focused on the rare cases where it fails. An im- 
portant example is the Kosterlitz-Thouless transition, 5 where 
no symmetry is broken and the phases are instead distin- 
guished by the asymptotic behavior of correlation functions. 
More recently, a class of quantum phase transitions has been 
proposed that have, by contrast, order in both phases with dis- 
tinct order parameters 

Another family of unconventional phase transitions, shar- 
ing features with both of these precedents, is now understood 
to exist in certain frustrated systems. Examples include the 
classical dimer model on the cubic lattice^ 5 -! and Neel order- 
ing of the Heisenberg model on the pyrochlore latticeP^Mem- 
bers of this family have also been pr edicte d in models of the 
magnetic materials known as spin icej 17 l 18 land there is numer- 
ical evidence of continuous transitions in a uniform applied 
magnetic field along a (100) crystal directiorp^" 2 ^ and in a 
nonuniform field with a helical structure in real spaceP 3 J 2 ^The 
common feature of these transitions is that they occur with the 
system subject to a constraint that prevents the occurrence of a 
simple thermally disordered state. The transition instead sep- 
arates a conventional (possibly ordered) phase below the crit- 
ical temperature from a so-called Coulomb phasepla species 
of classical spin liquid, 26 above. 

The Coulomb phase occurs when the low-energy configu- 
rations obey a topological constraint expressible as a lattice 



Gauss law. If states within the constrained sector are degen- 
erate, with all others gapped by at least A, then, for T <k A, 
the effective coarse-grained description is a noncompact U(l) 
gauge theory. A characteristic feature of the Coulomb phase 
is the presence of long-range dipolar correlations, which re- 
sult from long-wavelength fluctuations of the emergent gauge 
fieldP 

In physical systems, such an "accidental degeneracy" must 
be split at some smaller energy scale V <s A, either intrin- 
sically (e.g., due to further-neighbor interactions or quantum 
effects) or as the result of an external perturbation (e.g., ap- 
plied magnetic field). While the Coulomb phase remains for 
T » V, the fluctuations are suppressed for T <k V and the 
unique ground state, or one of a set of symmetry-related states, 
is selected. Even in the former case, when no symmetry is 
broken, the qualitative distinction between the two extremes 
of T/V imply that they must be separated by a phase transi- 
tion, at a temperature Tq ~ V, when the dipolar correlations 
of the Coulomb phase are lost. 

The correct description of the critical behavior at this tran- 
sition must include the long-wavelength gauge-field fluctu- 
ations of the Coulomb phase and their suppression at the 
transition, and therefore goes beyond the standard Landau 
paradigm. Critical theories for transitions in this family 
have bee n found u sing mappi ngs to conventional ordering 
transition ;] 1 m9 l 2 °l 27 -1 and analyses^ 2 ! 13 ! 23 * in terms of the Higgs 
mechanism! 2 ^ Any symmetry-breaking order that may ap- 
pear below Tq is more usefully viewed as a secondary 
phenomenon.^ (In fact, order and spin-liquid behavior are not 
mutually exclusive ! 26 l 29 j) 

Strongly correlated phases often host fractionalized excita- 
tions, such as spinons in frustrated quantum magnets ^ and 
fractional-charge quasiparticles in quantum Hall statesPsED 
In the Coulomb phase, they take the form of defects in the 
Gauss law constraint, or "charges", and occur with energy 
cost A. The topological nature of the constraint implies that 
they are conserved by coarse-graining, and correspond to 
monopoles in the effective gauge theory. In the case of spin 
ice, these monopoles (equivalent to spinons in a pseudospin 
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FIG. 1. Schematic phase diagram for a system with a continuous 
transition out of the Coulomb phase, as a function of temperature T 
and monopole fugacity z = e~ A/r , where A is the energy cost for a 
single monopole. The Coulomb phase is qualitatively distinct from 
the paramagnet only at z = (blue line), so the confinement tran- 
sition is an isolated point. It nonetheless influences properties in 
a broad region of the phase diagram and leads to universal scaling 
forms in its vicinity. The dotted line (T > T c ) indicates a crossover 
from Coulomb-like behavior to a conventional paramagnet, while the 
dash-dotted line (T < T c ) is either a transition (in a conventional 
universality class) or a crossover. Both have the form z ~ \T - Tel*, 
where is a crossover exponent. The arrow shows an example path 
as T is reduced at fixed A. 



description 32 33 ) carry not only fictitious gauge charge but also 
physical magnetic charge!^ 

Instead of characterizing the transition through the dipo- 
lar correlations, one can define an "order parameter" using a 
test pair of monopoles, i.e., a pair imposed on a state that is 
otherwise free of defects (see Section [II A\ . In the Coulomb 
phase at T > Tc, a test pair of monopoles with opposite gauge 
charge feels an entropic interaction oc \/R, at large separa- 
tion R (in three spatial dimensions, 3D). A finite free energy 
is required to separate them to infinite distance, and they are 
therefore deconfined. In the "conventional" phase at T < Tc, 
the fluctuations of the gauge field are suppressed, and the con- 
straint implies that separating defects leaves behind a trail of 
disturbance (see Fig. [2] for an illustration). This costs free 
energy proportional to its length, and hence causes an un- 
bounded potential oc which confines the defects. 

While a test pair of monopoles provides a criterion for the 
phases of the defect-free system, the distinction is lost when 
thermally excited monopoles are present. Such defects can 
"neutralize" the test monopoles and hence prevent the un- 
bounded interaction required for confinement. 35 The confine- 
ment transition therefore exists strictly in the constrained limit 
T/A — » 0, as illustrated in the schematic phase diagram of 
Fig. [T] and is replaced at nonzero defect density by either a 
crossover or a conventional transition. 

The present work shows that the unconventional critical 
point nonetheless has important consequences for behavior 
in the physical regime where defects are merely energetically 
suppressed, rather than absolutely forbidden. If monopoles 



cost energy A, so that their fugacity is z — e 



-A/T 



they will 



the critical point at T = Tc and z — 0, giving clear signatures 
of unconventional criticality even in phenomena at z > 0. 

We consider models with discrete defects, and focus on 
those where the degrees of freedom themselves are discrete 
(e.g., spin ice). The theory can also be applied in cases 
with continuous degrees of freedom but discrete defects car- 
rying nontrivial topological index (see Section |V Q . When 
the defect charge is continuous, the transitions can still be 
characterized through confine ment, but have more conven- 
tional Landau-type descriptions JlstJQI (As noted by Isakov et 
al.J^ continuous monopole charge implies that the correlation 
length diverges only algebraically with T/A in the Coulomb 
phase.) 

This work treats explicitly only classical statistical models, 
but the scaling theory is also valid for thermal phase transi- 
tions in quantum models. A similar analysis may prove ap- 
plicable to quantum phase transitio ns, s uch as those predicted 
in the quantum spin ice materials! 29 ! 42 ^ This work will fo- 
cus on 3D systems, although many of the present conclusions 
also apply in 2D, where other methods have been successfully 
applied.™ 

A brief account of the theoretical analysis and some of the 
numerical results have been presented previously. 24 Related 
work includes that of Castelnovo et al.p^who provide a qual- 
itative discussion of the crossovers between various regimes 
for a broader class of strongly constrained systems, as well 
as numerical results for certain 2D models. Bergman et alP^ 
predict a transition in a dimer model (closely related to one 
subsequently observed in the cubic dimer model 9 ), and briefly 
discuss the effects of nonzero monopole fugacity. 



Outline 



Scaling theory is applied to general confinement transitions, 
incorporating the effects of monopoles, in Section [IT] A mi- 
croscopic model describing the transitions of interest is in- 



troduced in Section II A and a mapping is presented in Sec- 
tion |IIB| that allows the standard results of scaling theory to 
be applied. The consequences for behavior at the critical point 



and within the Coulomb phase are detailed in Sections II C 
andUlD] 

Specific examples are provided by applications to two tran- 
sitions occurring within a model of spin ice: In Section[IlIj the 
Kasteleyn transition occurring in the presence of a uniform 
applied field is studied. This is at its upper critical dimen- 
sion and so exhibits logarithmic corrections to scaling, which 
are calculated and verified using numerical simulations. Sec- 
tion IV deals with a transition in the presence of a nonuniform 



occur with a nonzero densit y 37 " 39 ! for any z > 0. Scaling the- 
ory constrains physical quantities in the region surrounding 



applied field with a helical structure in real space. The univer- 
sality class is predicted to be that of the 3D XY model, and the 
critical exponents are shown to agree with established values 
for this class. 

A discussion of potential tests in experiment and in simula- 
tions of other models is given in Section[V] 
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II. GENERAL SCALING THEORY 

The behavior at nonzero monopole fugacity near confine- 
ment transitions can be understood using the general frame- 
work of scaling theory and the renormalization group (RG). 
Corresponding to any continuous phase transition is a fixed 
point of the RG transformation, and the critical properties are 
governed by the structure of the RG flow in its neighborhood. 3 

In the case of a confinement transition at zero monopole 
density, the topological constraint is preserved by an appro- 
priately chosen coarse-graining procedure, and so is inherited 
by the fixed-point theory. The behavior for small z near the 
critical point is therefore determined by the most relevant per- 
turbation that breaks the constraint. 

In other words, there exist one or more eigenoperators of 
the RG that are topologically forbidden when monopoles are 
absent. We seek the one with largest eigenvalue y z , defined 
such that its conjugate field z renormalizes to b y, z under a 
rescaling by factor b. The appropriate scaling field z is in fact 
given by exactly the monopole fugacity z, as demonstrated 



in Section II B The argument is quite general, with conse- 
quences that apply to a number of frustrated and constrained 
systems. 

If y z > 0, the perturbation of allowing monopoles is rele- 
vant, as will be the case in the examples treated in detail here. 
Since monopoles are a relevant perturbation at the Gaussian 
fixed point corresponding to the Coulomb phase (see Sec- 
tion II D i, irrelevance of monopoles at the critical point re- 



quires the presence of a phase transition at z > 0. The analysis 
to follow assumes that nonzero monopole fugacity is indeed 
relevant. (An alternative scenario is discussed briefly in Sec- 
tion[VB]) 



Coloring modelsS^ESl can be included if B( has multiple com- 
ponents, and so hosts monopoles of multiple flavors. 

The criterion for confinement involves inserting a pair of 
monopoles into an ensemble that is otherwise monopole-free. 
The corresponding partition function is 



E 



(2) 



where the sum is over configurations with monopoles of 
opposite charge at sites i and j, obeying div,< B = Sn 
The normalized partition function, 



St}- 



G m (r u ) 



z 



(3) 



where ry is the separation, describes the spatial distribution 
of test monopoles inserted into the constrained system. One 
can alternatively interpret U m (r) — - In G m (r) as the effective 
entropic interaction between the test pair, resulting from fluc- 
tuations of the gauge field. 

According to the definition of Section[TJ monopoles are "de- 
confined" if G m (r) has a nonzero limit as \r\ — > oo. In this 
case, U m (r) approaches a finite constant (with Coulomb-law 
corrections), and a finite energy is required to separate the 
monopole-antimonopole pair. In a confining phase, the po- 
tential instead grows without limit, and so G m (r) decreases 
exponentially to zero as \r\ — > oo. 



B. Monopole scaling field 



A. Microscopic model 

For the sake of concreteness, we focus on a discrete clas- 
sical model, although much of the analysis will apply to a 
broader class of systems. Suppose the system comprises de- 
grees of freedom Be defined on the links t of a lattice, scaled 
such that the monopole charge div, B takes integer values, 
where div, denotes the lattice divergence 49 at site i. The parti- 
tion function can then be written as 

z = JV-«>V^], (i) 

\B[) 

where a term depending on the monopole number has been 
separated from the rest of the action (i.e., configuration en- 
ergy) S[B]. Such a separation is unambiguous when the 
monopole energy cost is sufficiently large that the gap A above 
the constrained (zero-monopole) sector dominates other en- 
ergy scales, such as T and the splitting within this sector. Cor- 
rections to the harmonic form of the monopole energy, omit- 
ted from Eq. ([TJ, are negligible when z <K 1 and monopoles 
are rare [see Eq. ([6j, below]. 

Particular instances of Eq. ([TJ include the model of spin ice 
treated in Sections[TIl]and[lV] as well as other ice and current- 
loop models, as well as dimer models on bipartite latticesP^l 



We now turn to the problem of determining how the scaling 
field z is related to the microscopic parameter z- One requires, 
as usual, only the leading-order behavior close to the transi- 
tion, and only up to a constant of proportionality. 

To clarify the problem, it is worth contrasting with the more 
usual case where a perturbation appears as an additional lo- 
cal term in the microscopic Hamiltonian. The correspond- 
ing scaling field can then usually be determined simply by 
considering the transformation of this perturbation under the 
symmetries, and identifying the most relevant compatible RG 
eigenoperator. 3 In the present case, increasing z from zero in- 
stead reduces the energy cost of a monopole, previously infi- 
nite, to oc |lnz|. 

In fact, it will be demonstrated that the appropriate scaling 
field is simply the monopole fugacity z. This can be estab- 
lished using a mapping to an ordering transition, which re- 
places the topological constraint by a symmetry and allows the 
monopole fugacity to appear explicitly as an additional term 
in an effective action. While the mapping was originally used 
to relate the in teger current-loop model to the Villain form of 
the XY modelp2E21 the argument relies only on the universal 
features of the Coulomb phase, and so applies to a broad fam- 
ily of confinement transitions in constrained systems. 

Starting from Eq. ([TJ, one can introduce an explicit sum 
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(4) 



over the monopole numbers by writing 

Z- (do^ 

U [BMni] 

= f D6 X 6 " S[BH Zf ^ Smdf ° 6 ' Z ' "' e ' ' ^ 

^ IBf) {«,) 

where the angular variables 0, e [-n, n) constrain n, = div,- B. 
In Eq. ((5), the lattice form of the divergence theorem has been 
used; grad f is the lattice gradient on link I. 

In the case z — 0, the only nonzero terms are those with n, = 
everywhere. The remaining sum over B[ can in principle be 
carried out, leaving an effective action >S eff [#] that depends 
on the original action S[B]. This can be calculated explicitly 



only in simple cases (see Section IV A for an example), but 
the form of the coupling between and B( in Eq. ([5]) implies 
that it always involves only functions of grad f 6. It therefore 
has a global "XY" [or U(l)] symmetry under uniform shifts 
of 6i. 

The presence of an XY symmetry when z = suggests the 
possibility of ordering of the angle variables, and of a conven- 
tional phase transition from disorder to order as the parame- 
ters in S ef£ are varied. In fact, this transition corresponds to 
exactly the confinement transitions that we wish to describe: 
By writing G m in a form analogous to Eq. Q, one sees that it 
equals the (disconnected) correlation function (e~ ie, 'e +lfl (see 



Section II C 2 1. Deconfinement of test monopoles, signaled by 
a nonzero large-separation limit of G m (r), therefore coincides 
with long-range order in the angular variables. This identifies 
the Coulomb phase with the XY-ordered phase and hence the 
deconfinement transition with a conventional ordering transi- 
tion. (Note that the phases are "inverted" by the mapping, 
in the sense that the higher-temperature Coulomb phase maps 
to the ordered, and so lower-temperature, phase of the angle 
variables.) 

For < z <K 1, the sum over n, in Eq. |5]l can be evaluated 
separately at each site, giving 



JVe-"^ = l+ z (e- ie '-+e ie '-) + 



(6) 



cosOj 



The additional contribution to S eS of -2z2, cos ^i corre- 
sponds to an applied field /ixy k Z acting on the angle vari- 
ables. Increasing z from zero therefore explicitly breaks the 
XY symmetry, eliminating the possibility of an ordered phase. 
In these terms, monopole fugacity appears as the coefficient of 
an additional term in the action, and the XY symmetry allows 
the appropriate scaling field to be identified as simply /i X y ~ Z- 
It is important to note that neither the relationship between 
the XY-ordered and Coulomb phases nor the identification of 
the scaling field z relies on the details of the phase transition. 
There is no assumption, in particular, that the transition be- 
longs in the XY universality class. Indeed, when S[B] in- 
cludes interactions between B( on different links, the sum over 
Be will not factorize; the effective action .S eff [£>] may then in- 
volve complicated long-range couplings whose form will de- 
termine the universality class. 54 (A simpler case where the 



class is not XY is treated in Section III ) It is the generic nature 
of the Coulomb phase, and the relation between its topolog- 
ical order and the breaking of XY symmetry, that allows the 
scaling field to be identified. 



C. Scaling forms 

On the assumption that the monopole fugacity z is relevant 
(i.e., that y z > 0), the behavior in the neighborhood of the criti- 
cal point T — Tc and z = is described by the standard theory 
of crossover scaling with two relevant variables!^ This leads 
to scaling forms for thermodynamic quantities and correlation 
functions, expressed in terms of a small number of universal 
functions and critical exponents, determined by the properties 
of the critical fixed point. 

Consider, for example, the reduced (i.e., divided by T) free- 
energy density, defined by / = -A^ 1 lnJ2, where N is the 
number of degrees of freedom. The singular part f s of this 
function obeys 



Mt,z) 



j2-a 



® ± (z/\tf), 



(7) 



for sufficiently small t and z, where t = (T — Tq)/Tc is the 
reduced temperature. The critical exponents a and <p and the 
function cD ± are universal, in the sense that they depend only 
on the universality class of the transition. The subscript on i> ± 
indicates that the function may also depend on the sign of t, 
taking different forms above and below the critical tempera- 
ture. 

The scaling form in Eq. (j7]i follows directly from the fact 
that, in the neighborhood of the critical fixed point, rescaling 
by a factor of b replaces the values of t and z by tb y ' and zb y - 
respectively. 3 The values of the critical exponents are related 
to the RG eigenvalues y, and y z by 



d 

a = 2-~ 

yt 

yt 



(8) 
(9) 



Here a is the standard "specific heat" exponent, while <f> de- 
scribes the effect of nonzero monopole fugacity on the critical 
behavior. It will be referred to as a "crossover" exponent, 3 
since, as illustrated in Fig. [T] it controls how the system 
crosses over between different regimes when t and z are var- 
ied. In particular, nonzero z has greatest effect when \tf < z; 
otherwise the argument of ± is small and / s may be approx- 
imated by its z — behavior. 

It should be noted the expression Eq. |7]) for / s will in some 
cases need to be generalized to include more than two relevant 
fields. In particular, in cases where a symmetry is broken at 
the confinement transition, the order parameter couples to an 
additional scaling field, introducing (in most instances) a third 
independent critical exponent. Note also that scaling of f s 
applies even when there is a phase transition for z > 0, and 
that it then constrains the phase boundary 7c(z), as illustrated 
in Fig. [T] 

Scaling forms for other thermodynamic quantities, such as 
heat capacity, follow in the standard way from Eq. |7]) or 
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its generalizations, 
monopoles, 



Using Eq. (HI), the absolute density of 



J]\div i B\ 



(10) 



The relative sizes of these length scales also determine the 
crossover of the thermodynamic functions. The scaling form 
for f s in Eq. (|7]), for example, can be rewritten as a function 
of £/A m , so different regimes correspond to different relative 
magnitudes of the two length scales. 



can be expressed as a derivative of the (reduced) free-energy 
density, 



2. Monopole distribution function 



f, 



(11) 



when z is small (so that fluctuations with |div, B\ > 1 are rare). 
It therefore has a scaling form following from Eq. (|7J, 



p m (t,z)~\t\ 2 - a ®™(z/\t\*), 



(12) 



where (D™(x) = x<&' ± (x). (The monopole density, including 
any nonsingular part, vanishes at z = 0.) 

When testing the predictions of scaling theory, especially in 
numerical simulations, it is important to bear in mind finite- 
size effects. For a system of linear dimension L cc N l ^ d , in- 
tensive quantities such as the free-energy density will also be 
functions of L\t\ v , where v is the correlation-length critical ex- 
ponent, 

(13) 

y t a 



1. Correlation functions and length scales 

Similar scaling forms can be written for correlation func- 
tions. A two-point correlation function depends on the mag- 
nitude of the separation r, and may also have nontrivial sub- 
lattice and direction dependence. While scaling theory itself 
has little bearing on the latter, the effective theory describing 
a fixed point will often have sufficient symmetry that the sub- 
lattice and direction dependence can be strongly constrained. 

Scaling theory implies that, at large separation r — \r\, two- 
point correlation functions can be expressed in terms of com- 
binations r|f| v , rz v ^, and r/L. One is therefore led to define 
length scales £ ~ \t\~ l ^ y ' and A m ~ z -1 ^* associated with the 
relevant scaling fields t and z. The former takes the place of 
the correlation length 51 (which, according to a naive defini- 
tion, diverges throughout the Coulomb phase), while the latter 
gives the characteristic separation between monopole defects. 
Both £ and A m are infinite at the confinement-transition critical 
point. 

These scales govern crossovers between forms of the cor- 
relation function characteristic of different fixed points. For 
example, for t > and sufficiently small z (and L = oo), one 
can have a <k £ <k A m (where a is the lattice spacing). There 
will then be three distinct regimes of the correlation function: 
For a <K r <K one has critical correlations, governed by 
the critical point. For £ <K r <K A m , the correlations take 
the dipolar form characteristic of the Coulomb phase. Finally, 
for r » A m , one has conventional (exponential) paramagnetic 
correlations. 



While not expressible through correlations of local quanti- 
ties, the monopole distribution function G m , defined in Eq. ([3J, 
behaves analogously to the correlation functions discussed 
above. In fact, as noted in Section II B the mapping to an- 
gular variables gives the relation 



<e -is, e +ie y); 



(14) 



which expresses the distribution function as a standard cor- 
relation, and leads immediately to its scaling behavior. [Ex- 
pressing Eq. |2} in the same form as Eq. Q, Eq. ( 14 1 follows 
immediately. One can alternatively replace the fugacity by a 
nonuniform field z;, and take derivatives with respect to z; and 



As discussed in Section II A the transition separates two 



different asymptotic behaviors of G m (r), with a nonzero limit 
for \r\ — > in the Coulomb phase and exponential decay in the 
ordered phase. In the neighborhood of the transition, scaling 
theory giveP 



G m (r,t,z) ~ r 



,rz v/0 ) 



(15) 



where r m± is a universal function. (Possible anisotropy and 
sublattice dependence have been omitted.) 

This scaling form implies, in particular, that the monopole 
distribution function becomes a power law at the critical point, 
with exponent 2(d - y z ) = 2(d - <p/v). As demonstrated be- 



low in Section IV B this provides a method of determining 
the exponent (p in numerical simulations carried out at zero 
monopole fugacity. Such a calculation then provides quanti- 
tative predictions for simulations at z > 0, even for transitions 
such as in the cubic dimer model (see Section [VB] i, where the 
critical exponent <p is not known a priori. 



D. Scaling theory in the Coulomb phase 

The scaling theory described here can also be applied in the 
Coulomb phase itself, where the action density is given by the 
pure continuum U(l) gauge theory] 41 1 55 ! 56 1 



Xcouiomb = \k\B\ 2 = ^K\VxA\ 2 , 



(16) 



with K a positive constant. At this Gaussian fixed point, 
power counting shows that all analytic perturbations, corre- 
sponding to smooth fluctuations in the field B, are irrelevant. 
Monopole fugacity is relevant, however, and one can infer its 
(exact) RG eigenvalue = d from the nonzero limit of the 
test-monopole-pair correlation function G m in the Coulomb 
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phase. 57 Scaling, with the single scaling field z, then implies interaction J > 



z- l/d and 
Z. Taking the 



that the monopole separation obeys A„ 
the free-energy density scales as f s {z) ~ Z 
derivative with respect to lnz gives p m ~ z, as simple con- 
siderations would suggest. (Different forms hold in confined 
phases; see, for example, Section [Til C| ) 

Starting from the scaling form near the critical point, 
Eq. ( 12 1, and taking z — > at fixed t > 0, one should recover 
the Coulomb-phase result of monopole density proportional 
to z. This implies that the universal function O™ has limiting 
behavior <E>™(x) ~ x (for t > 0), and hence that the "monopole 
susceptibility" vanishes as dp m /dz\ : =o ~ t 1 ^"^ = f v(flL, -> when 
approaching the confinement transition. 

The fixed-point action Xcouiomb gives dipolar correlations of 
the field B, and hence "pinch points" in the structure factor for 
neutron scattering. 25 For z > 0, the asymptotic form of the cor- 
relations is an exponential decay with length scale A m ~ p„ , 
and the pinch points acquire a width ~ A^ . In the neighbor- 
hood of the critical point, this scales as /t m ' ~ z 1 ' f~^' , so the 
pinch points broaden as t — > + at fixed z > 0. 



III. KASTELEYN TRANSITION: SPIN ICE IN A (100) 
FIELD 



In this section and the next, the general scaling theory de- 
veloped in Section will be applied to two specific transi- 
tions occurring in a nearest-neighbor model of spin ice. The 
critical theory for both is well understood, allowing detailed 
predictions to be made for the behavior at nonzero monopole 
fugacity. 

This section treats the Kasteleyn transition that occurs in 
the presence of a uniform field along a (100) crystal direction. 
There will be considerable overlap between our results and 
those of previous workf^EH the emphasis here is on showing 
how they can be understood using scaling theory. 



A. Model: nearest-neighbor spin ice 



A realistic model of spin ice 17 18 involves unit-magnitude 
classical spins S( defined on the sites of a pyrochlore lattice, 
illustrated in Fig. |2] By identifying these with the links i of a 
diamond lattice, 25 one can cast this in the form of the general 
model Eq. (jlj, introduced in Section II A 

The spins are subject to a strong single-ion easy-axis 
anisotropy enforcing St - 2B(8(, where Be = ±i. The unit 
vector 6r points along the (111) direction of the diamond link 
£, with each link assigned a direction according to an arbi- 
trary convention. In the spin-ice materials such as Ho2Ti207, 
the anisotropy is effected by crystal fields that produce a split- 
ting above the Ising doublet of D > 100 KP-^ This is by far 
the largest scale in the problem, so the spins will be treated as 
binary degrees of freedom throughout. 

In a nearest-neighbor model of spin ice, a ferromagnetic 



Hi = 



(17) 



acts between neighboring pyrochlore sites (£€') or, equiva- 
lently, between all pairs of spins on each tetrahedron. With 
all the unit vectors 8( chosen to point towards a certain sublat- 
tice of the (bipartite) diamond lattice, 8c -dp = — I for every 
neighboring pair. The interaction Hj can therefore be rewrit- 
ten, by completing the square, as 



H^+y^PViB) 2 -!]. 



(18) 



The energy is minimized by any configuration of the spins 
where B{ is divergenceless^j on every diamond site. This 
constraint, more commonly expressed by saying that each py- 
rochlore tetrahedron has two spins po inting in and two point- 
ing out, is referred to as the "ice rule"p2EH 

The minimal defects in this constraint are diamond sites 
where div, B = ±1 (three spins out, one in; or vice versa), 
costing energy A = \j. (The link variables Be are scaled so 
that their lattice divergence div,- B takes integer values, as re- 



quired in Section II A ) These are monopoles in the sense that 
they carry "charge" in the lattice gauge theory. The Boltz- 
mann weight e~ Hj l T can be identified with the factor z^^'^ 
in Eq. |T|), so the monopole fugacity is z = e~ A/T = e - ?- 7 ' 7 . 

In the nearest-neighbor model, all states within the ice-rule 
manifold are exactly degenerate, so (assuming ergodicity) the 
system continues to fluctuate throughout this manifold and 
there is no ordering transition even in the limit T/J — > 0. 
When constrained in this way, the system exhibits a Coulomb 
phase: test monopoles are deconfined, with a Coulomb inter- 
action of entropic origin at large distance.^ 

In dipolar spin ice, realized in materials such as Ho2Ti2C>7, 
the spins are in fact subject to strong dipolar interactions, 
but the neglect of further-neighbor couplings is less severe 
than it may appearP^EH Within the ice-rule manifold, their 
only effect is a small splitting that causes (in simulations, 
at least) a first-order transition into an ordered state at very 
low temperature. 6 - 1 The magnetic dipole moment associated 
with the spins has the important effect that the monopoles of 
the effective gauge theory are in fact also physical magnetic 
monopoles, 34 and interact through a magnetostatic Coulomb 
force (besides the interactions induced by the fluctuating field 
B(). In this case, A includes a magnetostatic contribution from 



the field energy of an isolated monopole (see Section V A I 



B. Kasteleyn transition 

An applied magnetic field h couples to the spins through a 
Zeeman term, 



t 

= ~Yph-8 e )B e , 



(19) 
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FIG. 2. Left: Ground state of spin ice in a uniform field along the [001] crystal direction (vertical). Ising-like spins (green arrows) are arranged 
on the sites of a pyrochlore lattice, a network of comer-sharing tetrahedra. Each spin has a strong easy-axis anisotropy D forcing it to lie along 
the axis joining the centers of the two tetrahedra to which it belongs. The coupling / in Eq. selects a low-energy "ice-rule" manifold, 
where every tetrahedron has two spins pointing in and two pointing out. An external magnetic field h applied parallel to one of the cubic axes, 
with T <k \h\ <k J <k D, polarizes the spins, as shown. For T of order \h\, a Kasteleyn transition occurs, and strings of flipped spins proliferate, 
spanning the system in the direction of the field. Right: Strings terminate at defects in the ice rule, where a tetrahedron has three spins out and 
one in or vice versa [top (red) and bottom (blue) tetrahedra respectively]. The two phases of the system can be distinguished by the free-energy 
cost for separating a test pair of defects of opposite sign. In the saturated phase, this grows without limit as they are separated, so the defects 
are confined. 



and, for h = \h\ <K J, can be viewed as a perturbation split- 
ting the energies of the ice-rule states. For h small compared 
to T, the system continues to fluctuate within this manifold, 
but for sufficiently large h/T, the fluctuations are suppressed 
and a confinement transition occurs. In this section, we are 
concerned with a uniform field parallel to one of the cubic 
axes of the fee unit cell of pyrochlore, which is knowrPUl 
to cause a "Ka steleyn transition", with an unusual one-sided 
character! 45 ^ 61 1 

Consider a field h = V3/z u W[00i]> where K[ooi] is a unit vec- 
tor and the factor of V3 accounts for the cosine of the angle 
between the field and the local Ising axes. The ground state 
of H = Hj + Hh then has every spin aligned with K[ooi], to 
the extent permitted by the Ising anisotropy, as illustrated in 
Fig. [2] (left). The magnetization density, defined by 



42> 



(20) 



div B — 0, by flipping a closed loop of aligned spinsP^ In a 
fully polarized state, however, the only such loops span the 
whole system in the direction of polarization (assuming peri- 
odic boundaries in this direction). The minimal excitations 
above the ground state (in the limit J — oo) are therefore 
strings of flipped spins spanning the system in the [001] di- 
rection. 

Each flipped spin costs an energy 2\h ■ S(\ = -4= A, so the 
energy of a string is 2h n L^, in a system of linear extent (in 
the direction parallel to the field) L\\. Each string also con- 
tributes entropy L\\ In 2 (in the limit of low string density) due 
to the multitude of its possible paths. A single string therefore 
gives a contribution /\ tr i n g = L\\ (2h u - Tin 2) to the free en- 
ergy that is infinite in the thermodynamic limit L\\ — > oo and 
that changes sign at the temperature^ 



2 2 
In 2 V31n2 



h. 



(21) 



(where N is the number of spins), is saturated: m = (m) = 
m sat M[ooi] = 4|M[00i]- As the temperature is increased and the 
ratio h a /T is reduced, but remaining in the limit T <Z J, a 
transition takes place from this fully polarized state into the 
Coulomb phase. This transition, which breaks no symmetries 
and occurs when the system is strictly constrained to the ice- 
rule manifold, is closely analogous to the Kasteleyn transition 
occurring in dimer modelsP^lsD 

To understand its nature, consider the excitations above 
the ground state, which clearly involve flipped spins. To re- 
main within the ice-rule manifold, it is necessary to maintain 



For T < 7k, Turing is macroscopic and positive, so the num- 
ber of strings is precisely zero: all spins are aligned with the 
field and there are no fluctuations. Above 7k, the string den- 
sity increases continuously from zero, reducing the magneti- 
zation from its saturated value m sat and restoring the fluctu- 
ations that characterize the Coulomb phase. While no sym- 
metry is broken at the transition, an "order parameter" can be 
defined by taking the deviation from saturation magnetization, 
or equivalently the string (areal) density n = ^(m sat - m) = 
1(1 - m/m sat ), which vanishes in the low-temperature phase 
and increases continuously through the transition. 
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FIG. 3. Magnetization m versus reduced temperature t = (T-T K )/T K 
near the Kasteleyn transition at T = T K , calculated on the Bethe 
lattice. In the absence of monopoles (z = 0; top, black line), the 
magnetization is fixed at its saturation value m slit for all T < T K 
but decreases continuously for T > 7k. The transition is rounded 
for nonzero monopole fugacity z. For comparison, results of Monte 
Carlo simulations on the pyrochlore lattice (see Section [ill E 1[ > are 
shown with symbols; they exhibit the same qualitative features but 
agree quantitatively only far from the transition. (Finite-size effects 
reduce m from m sal at and just below 7k f° r Z = 0; see Fig. [5]) 



An alternative order parameter can be found by consider- 
ing a test monopole-antimonopole pair. Introducing such a 
pair into the saturated state requires flipping a string of spins 
joining the twoj^as illustrated in Fig.|5](right), and therefore 
costs free energy proportional to their separation. Monopoles 
are therefore confined in the low-temperature phase, and the 
transition to the Coulomb phase corresponds to deconfine- 
ment. This diagnostic for the transition agrees, as claimed in 
Section [I] with a definition based on thermodynamics or spin 
correlations. 

With a finite density of monopoles, i.e., for T/J > 0, the 
characteristic length for strings is set by the spacing A m be- 
tween monopoles rather than L\\, and F stl -i ns < oo even in the 
thermodynamic limit. The transition is then replaced by a 
crossover; at low temperature, a small density of short strings 
is thermally excited, reducing the magnetization from its sat- 
urated value. There is therefore an isolated critical point at 
T — T K at z — 0, with no phase transition at z > 0. 



C. Bethe lattice calculation 

The nearest-neighbor model of spin ice in a uniform field 
can in fact be solved exactly when the diamond lattice is re- 
placed by a Bethe lattice of the same coordination number. 
This amounts to a mean-field theory, and Monte Carlo sim- 
ulations on the original lattice show that it in fact provides a 
quantitatively accurate approximation for the latter, except in 
the neighborhood of transition, 19 as illustrated in Fig. [3] A 
brief description of the calculation is given in the Appendix [J 
see Ref.l2Tlfor details. 



Close to the critical point, the magnetization density can be 
expressed as 



m 



l-4z 



2/3 



"PI 



21n2 



-tz 



-2/3 



(22) 



where t — (T - 7\)/7\ is the reduced temperature, and ^(x) 
is defined as the positive solution of 



¥ 3 - x¥ - I = . 



(23) 



The singular part m s (t,z) = m(t,z) - m s . dt of the magnetiza- 
tion clearly obeys a scaling form similar to Eq. Q. The usual 
expression for the magnetization in terms of the free energy 
gives 



fs 



l,Z 



dt 



Muz) 



(24) 



t a <bf\zl\tf), 



where ( ± M) is a universal function, 
takes this form with ( ± m) (jc) = x^^i^x 
critical exponents <p = | and a — 0. (As remarked in Sec- 
tion IIC the subscript + provides a reminder the universal 



The Bethe-lattice result 

21 j^x' 2/i sgnf)] 2 and 



functions can also depend on the sign of t .) 

A quantity that will be particularly useful for the numeri- 
cal analysis of this transition (as well as for the helical-field 



transition of Section IV i is the variance of the magnetization, 
related to the "flux stiffness" in the Coulomb phase. 9 Specifi- 
cally, consider the component parallel to the field, defined by 



3N 
16" 



(25) 



which can be related to the susceptibility through the standard 
fluctuation-response relations, leading to 



Wn = 



1 / T 
41n2 



(26) 



The scaling dimension of W\\ is (within the Bethe-lattice cal- 
culation) given by —a = 0, which implies that the limiting 
form of W|| near the critical point is given by a universal func- 
tion of t/z 1 ^- This is confirmed in Fig. ffl which shows the 
full mean-field result for W\\ as a function of t and z- 

The absolute monopole density (per tetrahedron) can simi- 
larly be expressed as 



4z 4/3vj,j 



21n2 



-tz 



-2/3 



(27) 



which is of the form of Eq. ( 12 1, with the same exponent val- 
ues. Note that p m ~ zyft for z — > at fixed t > 0, in agreement 



with the general considerations of Section II D By contrast, 
p m ~ z 2 for t < 0; in the confined phase, monopoles are bound 
in charge-neutral pairs costing energy 2A and hence with ef- 
fective fugacity z 2 - 
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FIG. 4. Contour plot of uniform magnetization variance W\\, defined 
by Eq. l |25| l, calculated in mean-field theory (Bethe lattice). As the 
critical point (T = T K , z = 0) is approached, the contours take the 
form z~\T - r K |* (with cf> = |), as in Fig.[T| (Note that the contours 
are not evenly spaced.) 



These results provide an example of the general fact, proved 
in Section II B that z is the appropriate scaling field. As ex- 



pected for a mean-field theory, the critical exponents take ra- 



tional values. The quantum mapping of Section HID leads to 
a continuum critical theory for the Kasteleyn transition, which 
allows these values to be understood on the basis of Landau 
theory. 



D. Quantum mapping 

A continuum theory for the Kasteleyn transition can be 
found by exploiting a mapping to a model of quantum 
bosonsPSESThis identifies the strings of flipped spins, which 
traverse the 3D system in the direction of the field, with world 
lines of bosons in 2D space and imaginary time. The ab- 
sence of strings in the fully polarized phase maps to a vac- 
uum, while the higher-temperature Coulomb phase is a super- 
fluid of bosons, with a nonzero condensate order parameter. 
This mapping is clearly related to the general duality of See- 
but leads more naturally to a critical theory in the 



II B 



tion 

particular case of the Kasteleyn transition. 

More precisely, a transfer matrix T can be defined as the 
partial trace over all degrees of freedom lying between a pair 
of (001) planes, separated by distance <Sr. The full parti- 
tion function is X, = TrT L, ^ ST , and the thermodynamic limit 
L|| — > oo corresponds to the zero-temperature limit of the 
quantum problem. 4 One can then define the effective quantum 
Hamiltonian as ti — —St -1 InT. 

The precise form of *H is not important for the universal 
critical behavior, which can be inferred from the symmetries 
and topological properties of the two models. 1920 First, since 
strings span the system in the field (imaginary-time) direction, 
boson number is locally as well as globally conserved. The 
(2D) density of bosons is proportional to the string density n 
and hence to the deviation from saturation magnetization. In 



the original model, the applied field h u couples to the density 
of strings, so the chemical potential for bosons is // oc -h u . 

The Kasteleyn transition from the saturated to the Coulomb 
phase occurs when the applied field is reduced sufficiently that 
strings proliferate. (Since no closed loops are possible, strings 
either proliferate or are entirely absent.) This corresponds, 
in the quantum model, to increasing the chemical potential 
until the density of bosons becomes nonzero and they form 
a condensate. The quantum mapping implies that the critical 
behavior at this continuous transitio n 64 ! 65 ! is identical to that 
at the Kasteleyn transition. In these terms, the critical theory 
for the transition can be written down, using standard results, 4 
in terms of a critical field if/ representing fluctuations of the 
condensate order parameter. The action density isSHI 



£ K = ^5||iA+|V ± ^| 2 + ( 



+ !«!■ 



(28) 



where the coefficient (oc -/j) of the scalar \if/\ 2 tunes the system 
through the transition, so has been identified with t. Note the 
unusual sign of this coefficient: (if/) + gives the higher- 
temperature Coulomb phase. (The phases are inverted, as in 



the mapping of Section II B ) 

The term in Eq. (28 i involving a single derivative d\\ is al- 
lowed because reflection symmetry is broken by the presence 
of the uniform field, while a combination of reflection and 
spin inversion (time-reversal and particle-hole conjugation in 
the quantum problem) remains. 20 The presence of this sin- 
gle derivative implies that £ K has an invariance, at its criti- 
cal point, under an anisotropic scaling, where distances in the 
transverse directions are rescaled by b, while those parallel to 
the field are rescaled by b 2 . (In the quantum problem, 4 this is 
usually expressed by saying that the dynamical critical expo- 
nent is 3 = 2.) 

This has important consequences for the scaling theory of 
the Kasteleyn transition. First, it implies that correlations have 
anisotropic forms at the critical point. Second, results that in- 
volve hyperscaling, such as Eq. (|8j, are obeyed with effec- 
tive dimension d e s - 4, rather than the physical dimension 
d — 3. This further implies that the system is at its upper 
critical dimension,^ so one expects logarithmic corrections 3 
to results such as Eq. (|7]). 

Note that this is an example where the duality mapping of 



Section II B does not lead to a transition in the XY universality 
class. While condensation can be viewed as the ordering of 
angle variables (corresponding to the phase of if/), the resulting 
model has a continuum limit with single derivatives, which 
change the class of the transition. 



1. Effect of monopoles 

To determine the effect of nonzero monopole fugacity z on 
the effective quantum Hamiltonian fi, we return to its defi- 
nition in terms of the transfer matrix T. When z > 0, the 
partial trace involves configurations where a string terminates 
(see Fig. [2]), which occur with Boltzmann weight proportional 
to z, The matrix element between configurations on consec- 
utive (001) planes whose total string number differs by one, 
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denoted (rail \7~\ri), is therefore proportional to z (in the limit 
of small z). In terms of the effective quantum Hamiltonian, 
one has 



z oc (ra ± l|e _or/1 |ra> « (« ± 1|(1 - <Jr!W)|n> , (29) 

so (n ± 1 1*7/ |ra) ~ z. The contribution to the Hamiltonian then 
takes the form 



Hz = -J(.z)J](b l +bl), 



(30) 



where b t is the bosonic annihilation operator at site i on a (100) 
plane. This amounts to a field source J(z) ~ z, so the contri- 
bution to the continuum action is a term £ z z(i/r + 

This source term breaks boson-number conservation and 
hence the U(l) symmetry, corresponding to the XY symme- 



try appearing in the duality mapping of Section II B With 
this explicitly broken, there is no longer a phase transition 
for z > 0, in agreement with the microscopic considerations 
above. (The monopole-distribution function G m discussed in 



Section II C 2 can also be related to off-diagonal long-range 



order in the bosonic condensate. 1 ^!) 



E. Scaling theory and logarithmic corrections 

As the effective theory Xk + -Cz is at its upper critical di- 
mension, the Kasteleyn transition is governed by the Gaussian 
fixed point with t = z = u = 0. The RG eigenvalues for the 
scaling fields can be identified through dimensional analysis, 
which gives y, — 2 and y z = 3, along with y u = for the quar- 
tic coupling. Using the standard expressions for the critical 



exponents given in Section II C with an effective dimension- 
ality d e ff = 4, gives a — 2 and cf> — I, in agreement with the 



results of the Bethe-lattice calculation of Section III C Going 
beyond this mean-field theory, one expects logarithmic correc- 



tions to the scaling expressions of Section II C (as confirmed 
previously^ at z - 0), resulting from the marginally irrelevant 
coupling u. 

These corrections can be seen most clearly by focusing on a 
quantity with zero scaling dimension, which therefore remains 
finite at the transition within mean-field theory. A convenient 



choice is provided by W\\, defined in Eq. ( 25 1, which is propor- 
tional to the "flux stiffness", a quantity that has proven useful 



in analysis of related transitions (see also Section IV B I. Using 
Eq. ( |26| , it can also be related to the differential susceptibility 
(and hence to the heat capacity), which was used to demon- 
strate logarithmic corrections at z = in Ref.[T9l 

We follow the standard route to the leading logarithmic cor- 
rections to scaling, 3 which result from the slow decrease of 
the marginally irrelevant coupling u as the fixed point is ap- 
proached. One considers renormalization by a factor b » 1, 
chosen so that the characteristic length scale for fluctuations 
is reduced to the order of the lattice spacing. Following 
such a transformation, the effective value of u is given by 
Ub ~ (A In b) , where A is a positive constant independent (for 
sufficiently large b) of the original value of u. The parameters 
t and z and the field iff are similarly replaced by renormalized 



values tb ~ tb y ', Zb ~ zb y ', and \pb ~ ij/b** respectively, where 
the exponent x^ — 1 follows from dimensional analysis. 

At this point, the theory is sufficiently far from criticality 
that one can safely apply the results of Landau theory. These 
follow from the free-energy density, 



T K = -t 



+ l -u\ilf\ A - l -cz{^ + f), (31) 



where the derivatives in Eq. (28 1 have been dropped and c is a 



positive constant. The free energy is minimized by 

.1/3 / fM -l/3 \ 



(czY 13 / tu~ LID \ 

**+>(t.Z,H) = (-) ^(^J* (32) 



where *T is the function defined by Eq. ( |23) , 
One can write if/ ~ b~ x *if/ L . dn d RU (tb,Zb, u h ), and so 



czb 

Ub 



1/3 



tbU 



■1/3 



cAz In b 

zo 



1/3 



(czb) 2 ' 3 
(tA l '\\nby' 3 zf) 
fo(«F 3 



Using n ~ 4> and Eq. ( 26 1, 



xoz 



2/3 



(Inb) 



1/3 



(33) 
(34) 

(35) 



where ¥(x) - < F'(x) v F(x) and Wo and xq are unknown con- 
stants. 

In the present case, a suitable length scale is provided by 
b oc rT x l 2 , which gives the characteristic separation between 
strings. To arrive at analytically tractable expressions, but 
without affecting the result to leading logarithmic order, we 



replace n by the mean-field result n m f of Eq. ( 22 1, giving 

-l 



b = b 



,l/3xy ( 2in2 (z -2/3 



(36) 



where bo is a constant, used along with Wo and xo, to fit to 
numerical results. 

In the limit z — > 0, one finds (for t > 0) 



Wi ~ In 



to 



where to = |feyln2, consistent with the logarithmic singular- 
ity in the differential susceptibility observed at z = in earlier 
numerical simulations!^! 



1. Numerical results 

To test these predictions, we performed Monte Carlo (MC) 
simulations on the microscopic model Hj + Hh- In order to 
simulate the system efficiently in the regime of small or van- 
ishing z, a cluster (or "worm") algorithm was used, in which 
each MC step involves flipping a string of spins ! 66 ! 67 ! In the 
limit z — 0, only closed loops of spins are flipped, so the 
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FIG. 5. Variance W\\ of magnetization component parallel to applied 
field, illustrating scaling with logarithmic corrections near the Kaste- 
leyn transition. Symbols show the results of Monte Carlo simulations 
for a system of size L = 256 (with N = 16L 3 =s 4 x 10 6 spins), while 
the solid lines show a fit to Eq. {35} using scaling parameter b from 
Eq. {36}. There are significant deviations for the larger values of t 
and z, and also close to the transition for z = 0, where all strings 
span the system and so finite-size corrections become important. (In 
all other cases the characteristic monopole separation A m ~ z~ 1/3 is 
smaller than the system size, and so most strings are terminated by 
monopoles.) The dashed lines show the results of mean-field the- 
ory (Bethe lattice), which is quantitatively accurate only far from the 
critical point. The scaling dimension of W\\ vanishes, and so its mean- 
field value is finite as t — > + even for z = and L = <x>, in contrast 
to the logarithmic divergence of the full result. 



system remains within the ice-rule manifold. For z > 0, the 
algorithm also allows open strings, which can either trans- 
port monopoles or, with probability oc z 2 , create pairs of 
monopoles (of opposite sign). 

Previous work 19 ! 21 ! 22 ! using a similar algorithm has clearly 
demonstrated quantitative agreement between MC and Bethe- 
lattice results except quite close to the critical point, and log- 
arithmic corrections to scaling in its vicinity. Support for 



Eq. ( 35 I, which incorporates the effect of nonzero z near the 
critical point, is provided by Fig. [5] These results confirm 
the qualitative features predicted using the critical theory, viz. 
scaling with rational mean-field exponents up to multiplica- 
tive logarithmic corrections. Higher-order corrections are sup- 
pressed only by additional powers of logarithms and so are 
substantial, precluding accurate determinations of the param- 
eters or the critical exponents. 

The MC algorithm allows large system sizes to be simu- 
lated efficiently, particularly near the Kasteleyn transition, ob- 
viating the need for finite-size scaling: The logarithmic diver- 
gence of Wj| with t — > + occurs only in the thermodynamic 
limit when z — 0, and is cut off at a scale determined either by 
the system size L or by the monopole separation A m ~ z -1 ^ 3 . 
Even for the smallest exper imentally realistic values of z (of 



V A i, it is possible to simulate systems 



order 1CT 3 ; see Sectionl' 
with size L » A m . (In terms of the string picture described in 
Section |HI B| the absence of finite-size corrections for L » A m 
is understood by noting that nearly all strings are in this case 



terminated by monopoles, rather than spanning the system.) 
For the sizes used here, finite-size effects should therefore be 
significant only for z — 0, as indeed observed in the results 
shown in Fig. [5] 



IV. HELICAL-FIELD TRANSITION IN SPIN ICE 

As noted in Section [III] the nature of confinement transi- 
tions in spin ice due to an applied magnetic field depends on 
its orientation. A transition that is convenient from a theo- 
retical perspective, but considerably more challenging exper- 
imentally, can be induced by an applied field with a helical 
structure in real space! 23 ! 24 ! 

As in the Kasteleyn transition, the applied field selects a sin- 
gle configuration, so no symmetries are spontaneously broken 
in the low-temperature confining phase. In this case, fluctu- 
ations remain for all T > 0, in contrast to the fully saturated 
state occurring in the presence of a strong uniform field. In 
this regard, the helical-field transition is closer to a conven- 
tional ordering transition, and in fact belongs to the XY uni- 
versality class. A brief account of this analysis and some of 
the numerical results have been presented elsewhere! 2 ^ (Chen 
et alP^have studied a confinement transition in the same uni- 
versality class occurring in the "1GS" cubic dimer model.) 

We start once more from the nearest-neighbor model of spin 



ice introduced in Section III A with configuration energy H = 



Hj + Hi, as defined in Eqs. ( 17 1 and ( 19 1. Consider an applied 
field with uniform magnitude but a helical structure in space, 
as illustrated in Fig.pl The wavevector is q — — M[ooi]> such 
that the axis is aligned with the [001] cubic direction and the 
pitch is equal to the fee lattice constant a. The field is 



h(ri) = V3/! £ (cos q ■ n, sin q ■ rt, 0) , 



(37) 



where rt is the position of pyrochlore site I, with the origin 
chosen at the center of a tetrahedron. (The factor of V3 again 
accounts for the cosine of the angle between the local field 
and spin directions.) A transition occurs at a critical value of 
the ratio of the field strength h e to T, in the limit h e , T <K J. 

The presence of a transition, and the distinction from the 
Kasteleyn transition of Section[ni] can be appreciated by con- 
sidering the limit of a strong field and the unique state that 
results. The helical field in Eq. (37i is locally aligned with 
the [110] and [110] chains of the lattice and selects a spin 
spiral with the same structure, shown in Fig. [6] Such a con- 
figuration has zero uniform magnetization, and can therefore 
be connected to others within the ice-rule manifold by flip- 
ping spins along closed local (i.e., non-spanning) loops. Such 
loops cost Zeeman energy oc h e that is finite in the thermody- 
namic limit, and hence occur with nonzero probability for any 
T > 0. 

Above the transition, T > Tq, loops of flipped spins "pro- 
liferate", in the sense that a nonzero fraction span the system 
in the thermodynamic limit. 68 The effect of the transition on 
the confinement of monopoles can be understood, as in Sec- 
tion III B by noting that a test monopole-antimonopole pair 



must be joined by a string of flipped spins. This string grows 
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FIG. 6. Helical structure of the applied magnetic field and the 
unique spin configuration that minimizes the Zeeman energy. Spins 
are aligned along the [110] and [110] chains of the pyrochlore lat- 
tice (horizontal), following a right-handed helix with axis parallel to 
[001] (vertical). Unlike in the fully polarized state of Fig. [2] it is 
possible to flip a closed loop of spins and remain within the ice-rule 
manifold. The dashed line shows the shortest such loop, comprising 
8 spins. 



without limit as the pair are separated, costing free energy that 
remains finite only in the deconfined phase for T > Tq. 



A. Critical theory 

Confinement transitions within the zero-magnetization sec- 
tor of spin ice can be understood as Higgs condensa- 
tion of fictional "electric" charges, dual to the magnetic 
monopoles. 23 The helical-field transition naturally emerges as 
the simplest special case of this analysis, with a scalar (i.e., 
one-co mpone nt, complex) Higgs field. A standard duality 
mappingpEDthen implies that the transition belongs in the 3D 
XY universality class, as argued previously.^ 

In the particular case of the helical-field transition, a more 
direct route to a critical theory uses the mapping of Sec- 
tion |lIB]to express the spin-ice model in terms of angle vari- 
ables, as in Eq. d5]l. Since S[B] in this case involves a sum 
over terms acting on a single link, the effective action can be 
written as S eS = 2f -C(grad f 9) - 2z2/C° s #;- Proceeding to 
a continuum soft-spin description and expanding in powers of 
derivatives, the only quadratic derivative term that is consis- 
tent with the symmetrie s] 23 ! 70 ! is [V^| 2 , where iff ~ e ,e . 

The mapping of Section |HID| also applies in this case, with 
the resulting quantum model at half filling (one hard-core bo- 
son per two sites), since the uniform magnetization remains 
zero in both phases. The helical field maps to a staggered po- 
tential, and the ordered phase is a Mott insulator with the same 
symmetry as the potential's 3 ^ This transition is, as expected, 
in the 2 + ID XY universality classPl 



The fact that this transition belongs in a well-studied univer- 
sality class gives the advantage that precise values are known 
for the critical exponents. Using a combination of MC sim- 
ulations and series expansions, Campostrini et al. 73 found 
a = -0.0151(3) v = 0.6717(1), and = 0.3486(1), giving 
(f, = dv-/3= 1.6665(3). 



B. Numerical results 



The same MC algorithm described in Section III E 1 was 
used to study this model, in order to confirm the existence of 
the transition and its continuous nature, as well as its scaling 
properties both for z = and at z > 0. (Some of these numer- 
ical results have been presented in Ref. |24|) In this case, it is 
not possible to simulate such large sizes as for the Kasteleyn 
transition (presumably because fluctuations remain even in the 
low-temperature phase), and finite-size scaling is necessary to 
extract critical properties. 

As in the Kasteleyn transition, a useful quantity to locate 
and characterize the transition is provided by the variance of 
the uniform magnetization density (or, equivalently, the uni- 
form susceptibility). Because scaling is in this case isotropic, 
it is convenient to take the trace of the magnetization-variance 
tensor W(T, z, L) = 3L 4 (\m\ 2 ), where m is defined by Eq. (20 1 
and the choice of power of L will be explained below. In the 
constrained limit, z = 0, the magnetization variance is related 
to the "flux stiffness", which gives a measure of fluctuations 
between topological sectors. In the Coulomb phase, the effec- 
tive action of Eq. ( 16 1 implies 74 that W grows linearly with 
L, while the large fluctuations required to change topologi- 
cal sector are exponentially suppressed in the confined phase. 
One therefore expects a crossover asT/h s is increased through 
its critical value, which becomes sharper with increasing sys- 
tem size. 



Using the mapping described in Section II B one can in fact 
show that the uniform magnetization is conjugate to a twist in 
the boundary conditions of the angle variables {#,}. The scal- 
ing dimension of the latter must vanish, due to the (compact) 
U(l) symmetry of the confinement- transition fixed point. It 
follows that W(T,z,L) has zero scaling dimension (thanks to 
the appropriate power of L in its definition), 9 so its scaling 
form is given by 



W(T,z,L) ~ Q(fL 1/v ,zL wv ) 



r#/v\ 



(38) 



where Q is a universal function. This implies that plots of 
W(T, z = 0, L) versus T for different values of L should cross 
at 7c, as confirmed by the MC results shown in Fig. [7] Using 
the largest system sizes simulated, with L < 20 (i.e., 16L 3 < 
10 5 spins), we arrive at the value Tc/h e = 3.252(1) for the crit- 
ical temperature. The slope at the crossing, dW/dT\T=T c ,z=o, is 
furthermore predicted by Eq. (38 1 to be proportional to L 1/v , 
providing an estimate of v that is largely insensitive to the 
value of 7c. The power-law form is confirmed in Fig. |8j the 
fitted v is consistent with the 3D XY universality classP^I 



The scaling form of Eq. (38 1 is demonstrated for z > in 
Fig. [9] where W(T c ,z,L) and L- 1/v dW/dT\ T=Tc are shown to 
depend on z and L only through zZ/ /v . In both cases, convinc- 
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T/h e 



FIG. 7. Variance of uniform magnetization, W(T, z = 0, L), ver- 
sus T/h E for various system sizes L, showing a crossing at T c /h B = 
3.252(1), indicative of a continuous transition. The inset shows the 
region near the crossing, including larger system sizes (17 < L < 20) 
and error bars, both omitted from the main figure for clarity. The sys- 
tem comprises LxLxL fee unit cells, each containing 16 pyrochlore 
sites (spins). 
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FIG. 8. Determination of the correlation-length critical exponent v, 
using results at z = 0. Main figure: Log-log plot of dW(T, 0, L)/dT 
at T = Tc versus system size L, fit to oc L 1/r . The (blue) dashed line 
shows the best fit value of the exponent v = 0.670 ± 0.006, while the 
(black) solid and (purple) dash-dotted lines show fits with v fixed to 
its values for the 3D XY and Ising universality classes respectively. 
(The transition is predicted to be described by the former; the Ising 
class is displayed only for comparison.) The best fit v agrees with 
the XY class, but not Ising. Inset: Data collapse of W(T, 0, L) versus 
tL Uv using the XY exponent v = 0.6717 and T c /h s = 3.252. 



ing data collapse is found using the value <p — 1 .6665 of the 
3D XY universality class. Our best estimate of the exponent 
is (f> = 1.65(15), with a confidence interval based on the qual- 
ity of data collapse. Fig. 10 demonstrates the scaling form 

As noted in 
is linear for t > 0, demonstrat- 



given in Eq. ( fl"2| ) for the monopole density p„ 
Section III Dl the function O" 



ing deconfinement of monopoles in the Coulomb phase; it is 



quadratic in the confined phase, as in Section III C 

The critical exponent <p can also be determined using sim- 
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FIG. 9. Scaling at nonzero monopole fugacity z. Plot of W(Tc,z, L) 
(empty symbols, left scale) and 8W(T,z,L)/dT at T = T c (filled 
symbols, right scale) versus flf. The data lie on a single curve 
in each case, with v and (/> = 1.6665 taking values for the 3D XY 
universality class. (Colors indicate different values of L, using the 
same scale as the inset of Fig. [8]) 




FIG. 10 
fined in Eq 



Universal function O" 1 describing monopole density, de- 
The monopole density p m divided by \t\ 2 ~ a 



12 1. Ine monopole density p m divided by \t\- using 
the 3D XY exponent 73 a = -0.015, is plotted against z\t\~*. Data for 
various T and z (and fixed L = 16) collapse separately for each sign 
of t. (Symbols have the same meaning as in Fig.|9]) Inset: Same data 
on a logarithmic scale, showing that the universal function is linear 
for positive / (Coulomb phase, black symbols) and quadratic for neg- 
ative t (confined phase, red symbols). The upper (black) and lower 
(red) lines have slopes 1 and 2 respectively. 



ulations at z = 0, by measuring the test-monopole distribu- 
tion function G m , defined in Eq. Restricting to the critical 
point (t — and z — 0) and incorporating finite-size scaling 



replaces Eq. ( 15 i by 

G m (r,L) 



-2(d- 



y ' J T m (r/L) . 



(39) 



The partition function Xij in the presence of a monopole- 
antimonopole pair can be found in Monte Carlo simulations 
up to an L-dependent factor. We therefore calculate the ra- 
tio G m (R m - dx ,L)/G m (Rmia,L), with R max /L fixed and of order 
unity and R mm fixed and of order the lattice spacing. This ratio 
is proportional to LT^- d ~ y ^, allowing y z and hence <p = vy z to 
be found!^ 
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14<L<20 slope = -1.045 (6) 
12<L<20 -1.034(4) 



FIG. 11. Finite-size scaling of the test-monopole distribution func- 
tion G m (r, L) at the critical point, t = and z = 0. The value for 
maximal separation r = i? raax = V3aL/2 (with periodic boundary 
conditions) is normalized by its value for r = R m i n = = a/ V2 
(where e\ is one of the primitive unit vectors of the diamond lattice) 
and plotted as a function of system size L. Only even L can be used, 
to maintain a consistent definition of maximal separation, so smaller 
system sizes are included than in Fig. [8] The deviation of the point 
with L = 12 is likely due to finite-size corrections. Excluding it gives 
a best-fit line (solid) with slope -1.045(6), while including it (dashed 
line) gives a slope of -1.034(4). 



Using this procedure, illustrated in Fig. 11 we find 2(d - 
y z ) = 1.045(6). With v = 0.670(6) from Fig7|8] this gives the 
result <p = 1.660(15), consistent with the value <p = 1.6665(3) 
calculated using exponents reported in Ref.l73l 



V. DISCUSSION 

This work has presented a general theory to incorporate the 
effects of monopole defects on the critical behavior near con- 
finement transitions. The analysis applies to a broad family 
of phase transitions in frustrated systems that realize the so- 
called Coulomb phase, and leads to precise predictions that 
may be confirmed in numerical simulations or experiment. 

Two examples have been treated in detail here, both occur- 
ring in a model of spin ice. For the Kasteleyn transition in the 
presence of a uniform applied field, the logarithmic correc- 
tions to scaling have been calculated analytically and demon- 
strated using Monte Carlo simulations. A second transition, 
in the presence of a helical applied field, exhibits critical be- 
havior in the 3D XY universality class, as confirmed using 
simulations. 



A. Experiments: Spin ice 

The phenomenology of the Coulomb phase is well estab- 
lished in the classical spin-ice materials ) 17 ! 18 ! of which the 
most prominent examples are Ho2Ti207 and Dy 2 Ti207. In 
particular, heat-capacity and neutron-scattering experiments 



observe the residual entropy and dipolar correlations result- 
ing from a low-energy manifold that is extensively degener- 
ate but highly constrained. Experimental evidence also shows 
that the elementary excitations in this phase are single defects 
in the "ice rule" constraint (rather than single spin flips), and 
that these are in fact physical monopoles, carrying magnetic 
charge.^ 

No magnetic ordering is observed in these materials down 
to the lowest accessible temperatures (although a first-order 
transition is expected at a rather lower temperature due to 
long-range interactions between the magnetic moments^). In 
order to drive a transition, it is therefore necessary to apply 
external perturbations, such as pressure or, as in the models 



treated in Sections III and IV an applied magnetic field. The 



scaling theory presented here applies to any such transition, 
provided that it can occur within the constrained manifold and 
that it is continuous (rather than first order). 

Various measurements could provide evidence for the scal- 
ing forms presented in Section [II] including thermodynamic 
quantities such as heat capacity, magnetization, and uni- 
form susceptibility. If a phase transition remains at nonzero 
monopole density, scaling also constrains the form of the 
phase boundary as z — » 0, as illustrated in Fig. [T] Spin-spin 
correlations, measured in neutron scattering, 76 can also pro- 
vide evidence of scaling, as discussed in Section |IlC 1| 

In experiment, one may not have independent control over 
the perturbation V, such as in cases where it corresponds to a 
lattice distortion or further-neighbor interactions. One is then 
restricted to a one-dimensional path through the phase dia- 
gram of T/V versus z = e~ A / r , as illustrated in Fig.n] In other 
cases, such as when V is an applied magnetic field, the two- 
dimensional phase diagram can be explored by varying T and 
V separately. 

In the classical spin ice materials, dynamical freezing is 
observed on cooling to a temperature 7f ^ 0.6 K, causing a 
lower bound on accessible monopole fugacity of zt = e~ A/rf . 
Nonetheless, since strong evidence for the Coulomb phase 
exists in experiment, signatures of confinement transitions 
are likely also to be accessible. Indeed, taking the value 77 
A ^ 4K appropriate to Dy 2 Ti2C<7, one arrives at an estimate 
of Zf - 10~ 3 , indicating that the parameter regimes of Figs. [5] 
and|9]are within the reach of experiment. 

The case of most immediate experimental relevance is 
likely the Kasteleyn transition of spin ice in a (100) field, 
described in detail in Section III Experimental results for 
the ma gnetizat ion as a function of applied field have been 
reporte d 63 ^ 7 8 ! 79 ! an d are consistent with this and previous theo- 
retical workp2H3but confirmation of scaling as in Fig.[5]would 
require additional measurements in the temperature regime 
where the ice rules are quite well established, but the sys- 
tem remains ergodic. Logarithmic corrections are notoriously 
difficult to observe in experiment, but are possibly accessible 
in quantities such as the differential susceptibility or specific 
heat, which have zero scaling dimension. 

An important additional feature of the experimental sys- 
tems is the magnetic charge carried by monopoles,^ a con- 
sequence of the dipole moments of the spins. These cause 
Coulomb interactions between monopoles, in addition to 
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those induced (in the Coulomb phase) by fluctuations of the 
emergent gauge field. Since the scaling behavior at small z is 
governed by the response to an isolated monopole, the only 
effect is a finite magnetostatic contribution to A. It remains 
possible that the interactions decrease the window over which 
critical behavior is visible, as in the case of the liquid-gas crit- 
ical pointPMl] Such questions could be addressed by further 
numerical studies including Coulomb interactions. 

Realizations may also be provided by the rec ently discov- 
ered quantum analogs of the spin ice materials! 29 ! 42 ™ Con- 
tinuous thermal transitions in such systems are described by 
classical critical theories, as discussed here, and quantum fluc- 
tuations preserve dynamics even at very low densities of ther- 
mally excited monopoles. 



B. Cubic dimer model 



The classical cubic dimer model with interactions favor- 
ing columnar ordering provides an example of a symmetry- 
breaking transition from a Coulomb phase to a conventional 
ordered phase. 9 The balance of evidence supports the claim 
that the transition is continuous and presumably described 
by a noncompact Higgs model with SU(2) matter fieldsPSKl 
Testing the predictions of the present work in this context is 
likely to be more computationally demanding than in the tran- 
sitions studied here, but also more rewarding: It would pro- 
vide direct evidence for the unconventional nature of the tran- 
sition, and also demonstrate scaling in the presence of a tran- 
sition at z > (absent from both models simulated here). 

In this model, dimers are defined on the links of a cu- 
bic lattice, with a constraint demanding that each site be 
touched by precisely one dimer. On a bipartite lattice, the con- 
straint implies zero lattice divergence for an appropriately de- 
fined variable B(. &2 Defects in this constraint take the form of 
monomers, i.e., empty or multiply occupied sites, and increas- 
ing z from zero involves allowing these with nonzero Boltz- 
mann weight. 



As in Sections |III E 1 1 and [IV BJ one can characterize the 
transition by the flux stiffnessp^ defined in terms of the net 
"magnetization" of Bf . (For z > 0, this cannot be related to the 
number of dimers crossing a surface spanning the system, but 
is instead given by the staggered dimer occupation number.) 
This transition involves symmetry breaking and hence a con- 
ventional order parameter, whose Binder cumulant provides 
an additional means of characterizing the transition.™ When 
z > 0, both quantities are also functions of z/\tf. 

In contrast to the two examples studied here, this ordering 
transition is not in a universality class where the critical expo- 
nents are known accurately. A direct test of the scaling theory 
would be to find the exponent <p using the monomer-pair distri- 
bution function at z = 0, as in Section [TV B[ and then observe 
data collapse at z > using this exponent. 



Quadrupled monopoles: deconfined criticality 

The failure of the confinement criterion at z > and the 
relevance of z in the Coulomb phase do not in fact exclude 
the possibility that it is irrelevant at the confinement transi- 
tion. This scenario seems rather unlikely for single monopoles 



(especially given the duality mapping of Section II B i, but 
could occur in a model with a nonzero fugacity for multi- 
ple monopoles. For example, interpolating between square 
and triangular lattices by including dimers that span plaquette 
diagonals 83 84 is equivalent to allowing doubled monopoles. 

In fact, in certain quantum magnets,^ single and double 
monopoles (topological defects in the spin configuration, as 



discussed below in Section V C i are dynamically suppressed 
by Berry phases. Quadrupled monopoles are not suppressed, 
but are reckoned to be irrelevant, leading to so-called "decon- 
fined criticality". In the context of the constrained models that 
are the focus here, single and double monopoles can instead 
be explicitly forbidden (or directly suppressed by the micro- 
scopic Hamiltonian). Quadrupled monopoles can be included 
in the cubic dimer model through tetramers that occupy the 
four sites on a single cube sharing the same sublattice. 

When these defects are included, the Coulomb phase 
should be replaced by a (topologically ordered) paramag- 
net with exponential correlations, while the columnar-ordered 
phase remains. The phase transition between the two is, ac- 
cording to the deconfined criticality scenario, described by the 
same critical theory as the case without defects, viz. a non- 
compact Higgs model with SU(2) matter fields. A separate, 
more direct test of the irrelevance of quadrupled monopoles 
would involve measuring their distribution function, analo- 
gous to Eq. pi, as a function of the separation. 



C. Heisenberg model with suppression of "hedgehogs" 

Our results have been phrased in terms of discrete models 
with a local constraint, and the effect of a small density of de- 
fects in this constraint. Somewhat surprisingly, closely related 
physics can occur in an unfrustrated (classical) Heisenberg 
model, where the monopoles are instead topological defects 
in the spin configuration. When these discrete "hedgehog" 
defects are forbidden, 85 one finds a high-temperature phase 
that is described by the Coulomb phase of a noncompact U(l) 
gauge theory, and an unconventional transition into the mag- 
netically ordered phase. 

This description can be reached through the standard CP 1 
representation of the Heisenberg moments in terms of spinons 
minimally coupled to a compact U(l) gauge field. Hedgehogs 
correspond to monopoles, so their suppression is equivalent 
to the limit z — > of the discrete models treated here. The 
transition at which the Heisenberg spins order is described by 
(Higgs) condensation of spinons (electrically charged under 
the gauge field), and leads to confinement of test hedgehogs. 
Numerical results appear to confirm this picture of a continu- 
ous transition with unconventional critical behavior!^ 

Moving away from the limit where hedgehog configura- 
tions are completely forbidden amounts to increasing the 
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monopole fugacity z from zero. Scaling theory as presented in 
Section [II] should govern the behavior at z > near the order- 
ing transition, as can be tested through numerical simulations. 
In fact, because the critical theory is believed to be that of the 
Higgs transition of SU(2) matter fields, the universality class, 
and the exponent </>, should be the same as in the cubic dimer 
model. 



ACKNOWLEDGMENTS 

I am grateful to Claudio Castelnovo, John Chalker, Adrian 
Del Maestro, Michael Fisher, Ludovic Jaubert, Michael 
Levin, and Roderich Moessner for helpful discussions. This 
work was supported by JQI-NSF-PFC and AFOSR-MURI. 



Appendix: Bethe lattice calculation 

The nearest-neighbor model of spin ice can be solved ex- 
actly, including for z > 0, when the pyrochlore lattice is re- 
placed by a Husimi cactus constructed from tetrahedraP^ A 
brief outline of this solution is presented here; readers are re- 
ferred to Ref.|2T|for more details. 

The sites of pyrochlore are equivalent to the links of a dia- 
mond lattice, and a Husimi cactus of tetrahedra can similarly 
be constructed from a Bethe lattice with coordination 4. The 
absence of loops in the Bethe lattice allows the partition func- 
tion to be expressed exactly using a recurrence relation. We 
will consider a lattice of finite size with an open boundary, and 
subsequently take the thermodynamic limit of observables de- 
fined deep within the interior of the lattice. 

Consider first a single branch of the Bethe lattice, based 
at a "root vertex" with 3 neighbors. The fourth link of the 
root site has fixed spin cr, where <x = ±1 means the spin is 
aligned or antialigned with the applied field. Let Z^ a be the 
partition function for such a branch, with total depth N; the 
full lattice can be constructed as two such branches, of length 
N and N - I. The full partition function for a Bethe lattice of 
depth N is given by tracing over the common link, 



(A.l) 



and the magnetization (in the field direction) of the central 
spin can be written 



Mn — m S atZ N (Z^ + Zat_i i+ - Zw-Zjv_i _) ■ 



(A.2) 



By considering all configurations of the four links of the 
root vertex, one finds recursion relations 



Ziv+i,± - g Tl Zfl ± + 2Z\j^Z^ ± + 2g* 2 zZ^Z^ 



+ g^zZ 2 N± Z NT +g ± kz ] 



■jvt + z Z N ^Z N± , (A.3) 



where the Boltzmann weight g — e 2h ^ T = 2 T ^ T is split 
equally between the two vertices to which each spin belongs. 

There is no finite limit as N — > oo of Z^ a , but the ratio 



1/2 Zn- 



(A.4) 



(with the factor of g 1 included for convenience) has fixed 
points Y^ given by roots Y of the quartic equation 

zY 4 +(2+z 4 -g)Y 3 +3z(l-g)Y 2 +(l-2g-gz 4 )Y-zg = 0. (A.5) 
At this fixed point, the magnetization of the central spin, 



M = nu 



- Y 1 
+ Y 2 ' 



(A.6) 



is taken as representative of the bulk magnetization density m. 
For z = 0, there is a nontrivial solution Y = y[y, where 



2g-l 
2-z 



(A.7) 



as long as y > (g > I, T > T K ); otherwise the only solution 
is Y = 0. Near the critical point (at y = 0, z = 0), Y, z, and y 



are all small, so Eq. (A.5 i can be rewritten in terms of y and 
replaced by 



1 

yY--z 



0, 



(A.8) 



where terms of order y 5 ^ 2 and zy have been dropped. Expand- 
ing y - 2ii2f + 0(t 2 ) and using Eq. (A.6 1 gives the leading- 
order behavior expressed in Eq. (22 1. 



A similar enumeration of all configurations leads to 



Pm,N+l = 4Z Arllfe? ^Z^ + Za,_ + Zg +tl Z 3 N _Z N+ 

+ z A Z 2 N+ Z 2 N _) (A.9) 



for the mean absolute monopole charge at the central vertex of 



the Bethe lattice, defined as in Eq. ( 10 1. In the limit N — > oo, 
this becomes 



Az 



Y + z 3 Y 2 + Y 3 



1 + 4zY + 2(2 + z 4 )Y 2 + \zY 3 + Y 4 



(A. 10) 



which reduces to Eq. ( 27 1 near the transition 
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